home *** CD-ROM | disk | FTP | other *** search
/ TPUG - Toronto PET Users Group / TPUG Users Group CD / TPUG Users Group CD.iso / PET / S-Super PET / (s)t2.d64 / REGRESSION.FTN < prev    next >
Text File  |  2009-01-18  |  18KB  |  557 lines
  1. *          Program reg.test
  2. *  To test the polynomial fitting subroutine 'reg'.
  3. *
  4.       character fnamein , fnameout , pname
  5.       integer run
  6.       real x(15),y(15),a(11)
  7.       print,"Enter input data file name"
  8.       read,fnamein
  9.       print,"Output file name"
  10.       read,fnameout
  11.       open(unit=1,file=fnamein)
  12.  
  13. *     *** Read the run number and problem name
  14.  
  15.       read(1,1) run , pname
  16.     1 format(i2,a16)
  17.  
  18. *     *** Read the number of data points and degree of
  19. *     *** polynomial to be fitted.
  20.  
  21.       read(1,2) n , m
  22.     2 format(2i3)
  23.  
  24. *     *** Read list of dependent and independent variables.
  25.  
  26.       read(1,3) (x(i),y(i),i=1,n)
  27.     3 format(2f10.4)
  28.  
  29.       close(unit=1)
  30.       call reg(pname,n,m,run,x,y,a,2,fnameout)
  31.       end
  32.  
  33.  
  34.       subroutine reg (pr , n , m , nplot , x , y , a , lfn , fname)
  35. *
  36. *  Purpose - To perform polynomial regression
  37. *
  38. *  Usage -
  39. *            call reg (pr,pri,n,m,nplot,x,y,a,lfn,fname)
  40. *
  41. *  Arguments -
  42. *
  43. *      pr......Alphanumeric problem name.
  44. *      n.......Number of observations.
  45. *      m.......Highest degree polynomial specified.
  46. *              If best fit polynomial is of degree
  47. *              less than  m,  m  will be reset to
  48. *              degree of best fit polynomial.
  49. *      nplot...Option code for production of printer
  50. *              plots of regression results.
  51. *                  0  if plot not desired.
  52. *                  1  if plot desired.
  53. *      x.......vector of independent variables.
  54. *      y.......vector of dependent variables.
  55. *      a.......vector of regression coefficients.
  56. *              There will be  m + 1  entries in  a .
  57. *      lfn.....Logical unit number for output file.
  58. *      fname...Output file name
  59. *
  60. *  Subroutines required - pplot , orthls , coeffs , fity
  61. *
  62. *  Remarks -
  63. *      (1) The largest degree polynomial which may be fitted is 9.
  64. *      (2) The number of data points must be greater than the
  65. *          largest degree polynomial fitted.
  66. *      (3) A maximum of 50 observations is assumed.
  67. *
  68. *--------------------------------------------------------------------
  69. *
  70.       character pr , fname
  71.       real x(n) , y(n) , a(11)
  72.       real w(50) , t1(50) , t2(50) , t3(50) , plt (150) , z(50)
  73.       real c(10) , alpha(10) , beta(10) , b(10)
  74. *
  75. *--------------------------------------------------------------------
  76. *
  77.        open(unit=lfn,file=fname,access="sequential",recl=80)
  78. *
  79.      3 format("1Polynomial Regression.....",a16/)
  80.      4 format(" Number of Observations",i6//)
  81.      5 format(" Polynomial Regression of Degree",i3/1x,
  82. &              "***********************************"/)
  83.      6 format("   Coefficients"/3(1x,4e15.5/))
  84.      8 format("0","      --- Analysis of Variance ---"/)
  85.      9 format(" ",5x,"Source",8x,"DF",8x,"SS",11x,"MS",12x,"F"/)
  86.     10 format(" ",2x,"Regression...",i6,3f13.3)
  87.     11 format(3x,"Error........",i6,2f13.3)
  88.     12 format(3x,"Total........",i6,f13.3/)
  89.     13 format("   No Improvement")
  90.     14 format(" "//27x,"Table of Residuals"//" Obs. No.",5x,
  91. &             "X Value",6x,"Y Value",4x,"Y Estimate",4x,
  92. &             "Residual"/)
  93.     15 format(1x,i6,4f13.4)
  94.     16 format(" ",11x,"Improvement in SS",f16.3///)
  95. *
  96. *--------------------------------------------------------------------
  97. *
  98. *     *** Print problem name and no. of obs., check order
  99. *
  100.       write(lfn,3) pr
  101.       write(lfn,4) n
  102.       if (m .gt. 9) m = 9
  103. *
  104. *     *** Find mean and total sum of squares for dependent variable
  105. *
  106.       ybar = 0.0
  107.       ssat = 0.0
  108.       do 100 i = 1 , n
  109.          ybar = ybar + y(i)
  110.          ssat = ssat + y(i) * y(i)
  111.   100 continue
  112.       fn = n
  113.       ybar = ybar / fn
  114.       ssat = ssat - fn * ybar * ybar
  115.       nt = n - 1
  116. *
  117. *     *** Generate orthogonal polynomials and coefficients
  118. *
  119.       call orthls(x,y,w,n,0,0,c,alpha,beta,m,m+1,t1,t2,t3,nd1)
  120. *
  121. *     *** Check fit for successive degrees, and find sum of squares
  122. *         attributable to error.  Compute -
  123. *           ssar  - S of S attributable to regression
  124. *           smar  - Mean square for regression
  125. *           smae  - Mean square for error
  126. *           fval  - F-Ratio
  127. *           sumip - Improvement in SS
  128. *           nr    - D of F for regression
  129. *           ne    - D of F for error
  130. *           nt    - D of F Total
  131. *
  132.       sum = 0.0
  133.       do 200 k = 1 , m
  134.          k1 = k
  135.          call fity(x,n,0,c,alpha,beta,k,k+1,plt,t1,t2,iend3)
  136.          ssae = 0.0
  137.          do 120 i = 1 , n
  138.             temp = y(i) - plt(i)
  139.             ssae = ssae + temp * temp
  140.   120    continue
  141.          nr = k
  142.          ne = nt - nr
  143.          ssar = ssat - ssae
  144.          smar = ssar / float(nr)
  145.          smae = ssae / float(ne)
  146.          fval = smar / smae
  147.          sumip = ssar - sum
  148.          sum = ssar
  149.          write(lfn,5) k
  150.          if (sumip) 140 , 140 , 150
  151.   140       write(lfn,13)
  152.             go to 210
  153.   150    call coefs(0,c,alpha,beta,k,k+1,a,t1,t2,t3,ind2)
  154.          k2 = k + 1
  155.          write(lfn,6) (a(l),l=1,k2)
  156.          write(lfn,8)
  157.          write(lfn,9)
  158.          write(lfn,10) nr , ssar , smar , fval
  159.          write(lfn,11) ne , ssae , smae
  160.          write(lfn,12) nt , ssat
  161.          write(lfn,16) sumip
  162.   200 continue
  163. *
  164. *     *** Plot or not
  165. *
  166.   210 if (nplot) 400 , 400 , 220
  167. *
  168.   220 call fity(x,n,0,c,alpha,beta,k2-1,k2,z,t1,t2,iend3)
  169.       do 230 i = 1 , n
  170.          plt(i) = x(i)
  171.          i1 = n + i
  172.          plt(i1) = y(i)
  173.          i1 = i1 + n
  174.          plt(i1) = z(i)
  175.   230 continue
  176. *
  177. *     *** Print table of residuals
  178. *
  179.       write(lfn,3) pr
  180.       k = k2 - 1
  181.       write(lfn,5) k
  182.       write(lfn,14)
  183.       do 250 i = 1 , n
  184.          resid = y(i) - z(i)
  185.          write(lfn,15) i,x(i),y(i),z(i),resid
  186.   250 continue
  187.       call pplot(nplot,plt,n,3,3*n,n,1,68,lfn)
  188.   400 m = k
  189.       close(unit=lfn)
  190.       return
  191.       end
  192.  
  193.  
  194.       subroutine orthls (x,y,w,n,l,j,c,alpha,beta,k,k1,t1,t2,t3,ind1)
  195. *--------------------------------------------------------------------
  196. *     This subroutine computes the coefficients of the polynomial
  197. *     equation of degree k and the alpha and beta parameters.
  198. *--------------------------------------------------------------------
  199.       real x(n),y(n),w(n),c(k1),alpha(k1),beta(k1),t1(n),t2(n),t3(n)
  200. *--------------------------------------------------------------------
  201. *     Program initialization.
  202. *--------------------------------------------------------------------
  203.       kj1 = k - j + 1
  204.       if (kj1 .le. 0) go to 16
  205.       sum = 0.0
  206.       if (l .eq. 1) go to 3
  207.          do 2 i = 1 , n
  208.             t3(i) = x(i)
  209.             if (j .gt. 0) go to 1
  210.                sum = sum + 1.0
  211.                go to 2
  212.     1       sum = sum + x(i)**(2*j)
  213.     2    w(i) = 1.0
  214.          go to 7
  215.     3 do 6 i = 1 , n
  216.          t3(i) = x(i)
  217.          if (j .gt. 0) go to 4
  218.             sum = sum + w(i)
  219.             go to 5
  220.     4    sum = sum + w(i) * x(i)**(2*j)
  221.     5    x(i) = w(i) * x(i)
  222.          y(i) = w(i) * y(i)
  223.     6 continue
  224.     7 b = 0.0
  225.       r0 = sum
  226.       do 9 i = 1 , n
  227.          if (j .gt. 0) go to 8
  228.             t2(i) = 1.0
  229.             go to 9
  230.     8    t2(i) = t3(i)**j
  231.     9 t1(i) = 0.0
  232. *--------------------------------------------------------------------
  233. *     Begin computation.
  234. *--------------------------------------------------------------------
  235.       ii = 1
  236.    10 s = 0.0
  237.       do 11 i = 1 , n
  238.    11 s = s + y(i)*t2(i)
  239. *--------------------------------------------------------------------
  240. *     Computation of a coefficient in the polynomial equation.
  241. *--------------------------------------------------------------------
  242.       c(ii) = s / r0
  243.       if (ii .ge. kj1) go to 15
  244. *--------------------------------------------------------------------
  245. *     Computation of an alpha for the polynomial equation.
  246. *--------------------------------------------------------------------
  247.       sumxps = 0.0
  248.       do 12 i = 1 , n
  249.    12 sumxps = sumxps + x(i) * t2(i) * t2(i)
  250.       alpha(ii) = sumxps / r0
  251. *--------------------------------------------------------------------
  252. *     Computation of a new polynomial.
  253. *--------------------------------------------------------------------
  254.       do 13 i = 1 , n
  255.          temp  = t2(i)
  256.          t2(i) = (t3(i) - alpha(ii)) * t2(i) - b * t1(i)
  257.          t1(i) = temp
  258.    13 continue
  259. *--------------------------------------------------------------------
  260. *     Computation of a beta for the polynomial.
  261. *--------------------------------------------------------------------
  262.       r = 0.0
  263.       do 14 i = 1 , n
  264.    14 r = r + w(i) * t2(i) * t2(i)
  265.       beta(ii) = r / r0
  266.       r0 = r
  267.       b  = beta(ii)
  268.       ii = ii + 1
  269.       go to 10
  270. *--------------------------------------------------------------------
  271. *     Successful return.
  272. *--------------------------------------------------------------------
  273.    15 ind1 = 1
  274.       return
  275. *--------------------------------------------------------------------
  276. *     Error return.  Set all c coefficients, alpha and beta to zero.
  277. *--------------------------------------------------------------------
  278.    16 do 17 ii = 1 , k
  279.          c(ii) = 0.0
  280.          alpha(ii) = 0.0
  281.          beta(ii)  = 0.0
  282.    17 continue
  283.       c(k+1) = 0.0
  284.       ind1 = -1
  285.       return
  286.       end
  287.  
  288.  
  289.       subroutine coefs(j,c,alpha,beta,kc,k1,a,t1,t2,t3,ind2)
  290. *-------------------------------------------------------------------
  291. *     This subroutine computes the  a  coefficients for a
  292. *     polynomial of degree  kc  where  kc  is less than or
  293. *     equal fo  k .
  294. *-------------------------------------------------------------------
  295.       real c(k1),alpha(k1),beta(k1),a(k1),t1(k1),t2(k1),t3(k1)
  296. *-------------------------------------------------------------------
  297. *     Program initialization.
  298. *-------------------------------------------------------------------
  299.       kcj1 = kc - j + 1
  300.       if (kcj1 .le. 0) go to 9
  301.       b = 0.0
  302.       do 1 nn = 1 , kcj1
  303.          a(nn)  = c(nn)
  304.          t1(nn) = 0.0
  305.          t2(nn) = 0.0
  306.          t3(nn) = 0.0
  307.     1 continue
  308.       if (kc .le. j) go to 5
  309.       ii = 2
  310. *-------------------------------------------------------------------
  311. *     Begin computation.
  312. *-------------------------------------------------------------------
  313.     2 t2(ii) = 1.0
  314.       do 3 nn = 2 , ii
  315.          t3(nn)  = t2(nn-1) - t2(nn) * alpha(ii-1) - b * t1(nn)
  316.          a(nn-1) = a(nn-1) + c(ii) * t3(nn)
  317.     3 continue
  318.       if (ii .ge. kcj1) go to 5
  319. *-------------------------------------------------------------------
  320. *     Reset the vectors for the next coefficient.
  321. *-------------------------------------------------------------------
  322.       do 4 nn = 1 , ii
  323.          t1(nn) = t2(nn)
  324.          t2(nn) = t3(nn)
  325.     4 continue
  326.       b  = beta(ii-1)
  327.       ii = ii + 1
  328.       go to 2
  329.     5 if (j .le. 0) go to 8
  330. *-------------------------------------------------------------------
  331. *     Arrange coefficients properly if  j  is non zero.
  332. *-------------------------------------------------------------------
  333.       do 6 nn = 1 , kcj1
  334.          n1 = kcj1 - nn + 1
  335.          n2 = n1 + j
  336.          a(n2) = a(n1)
  337.     6 continue
  338.       do 7 nn = 1 , j
  339.     7 a(nn) = 0.0
  340. *-------------------------------------------------------------------
  341. *     Successful return
  342. *-------------------------------------------------------------------
  343.     8 ind2 = 2
  344.       return
  345. *-------------------------------------------------------------------
  346. *     Error return, set all  a  coefficients to zero
  347. *-------------------------------------------------------------------
  348.     9 do 10 nn = 1 , kc
  349.    10 a(nn) = 0.0
  350.       a(kc+1) = 0.0
  351.       ind2 = -2
  352.       return
  353.       end
  354.  
  355.  
  356.       subroutine fity(xf,m,j,c,alpha,beta,kf,kf1,yf,t1,t2,ind3)
  357. *---------------------------------------------------------------------
  358. *     This subroutine computes the fitted values for a
  359. *     given set of arguments with a polynomial of degree
  360. *     kf where kf is less than or equal to k.
  361. *---------------------------------------------------------------------
  362.       real xf(m),c(kf1),alpha(kf1),beta(kf1),yf(m),t1(m),t2(m)
  363. *---------------------------------------------------------------------
  364. *     Program initialization.
  365. *---------------------------------------------------------------------
  366.       kfj1 = kf - j + 1
  367.       if (kfj1 .le. 0) go to 7
  368.       b = 0.0
  369.       do 2 i = 1 , m
  370.          yf(i) = 0.0
  371.          if (j .gt. 0) go to 1
  372.             t2(i) = 1.0
  373.             go to 2
  374.     1    t2(i) = xf(i)**j
  375.     2 t1(i) = 0.0
  376. *---------------------------------------------------------------------
  377. *     Begin computation.
  378. *---------------------------------------------------------------------
  379.       ii = 1
  380.     3 do 4 i = 1 , m
  381.     4 yf(i) = yf(i) + c(ii)*t2(i)
  382.       if (ii .ge. kfj1) go to 6
  383. *---------------------------------------------------------------------
  384. *     Computation of new polynomial.
  385. *---------------------------------------------------------------------
  386.       do 5 i = 1 , m
  387.          temp = t2(i)
  388.          t2(i) = (xf(i) - alpha(ii)) * t2(i) - b * t1(i)
  389.          t1(i) = temp
  390.     5 continue
  391.       b  = beta(ii)
  392.       ii = ii + 1
  393.       go to 3
  394. *---------------------------------------------------------------------
  395. *     successful return
  396. *---------------------------------------------------------------------
  397.     6 ind3 = 3
  398.       return
  399. *---------------------------------------------------------------------
  400. *     Error return. Set all fitted values equal to zero.
  401. *---------------------------------------------------------------------
  402.     7 do 8 i = 1 , m
  403.     8 yf(i) = 0.0
  404.       ind3 = -3
  405.       return
  406.       end
  407.  
  408.  
  409.       subroutine pplot(n0,a,n,m,nm,nl,ns,nchar,lfn)
  410. *---------------------------------------------------------------------
  411. *       Purpose -
  412. *            Plot several cross-variables versus a base
  413. *            variable in strip chart form.  This is not
  414. *            a tidy program and may have bugs in the y-axis
  415. *            (looking sideways at the plot) labels.  Main use
  416. *            was for quick looks at long strip chart like
  417. *            arrays of data.  This program was evidently
  418. *            taken from an old IBM SSP package.
  419. *
  420. *       Usage -
  421. *            call pplot(n0,a,n,m,nm,nl,ns,nchar,lfn)
  422. *
  423. *       Arguments -
  424. *            n0      - Chart number (3 digits maximum)
  425. *            a       - Matrix of data to be plotted.  First
  426. *                      column represents base variables and
  427. *                      successive columns are the cross-variables
  428. *                      (maximum is 9)
  429. *            n       - Number of rows in matrix a
  430. *            m       - Number of columns in matrix a
  431. *                      (equal to the total number of variables)
  432. *                      Maximum is 10
  433. *            nm      - The total number of elements in a (usually n*m)
  434. *            nl      - Number of lines in the plot.  If 0 is
  435. *                      specified, 50 lines are used
  436. *            ns      - Code for sorting the base variable data
  437. *                      in ascending order
  438. *                        0  sorting is not necessary (already in
  439. *                           ascending order)
  440. *                        1  sorting is necessary
  441. *            nchar   - Number of print positions available
  442. *            lfn     - Logical file number for printed output
  443. *---------------------------------------------------------------------
  444.       real      ypr(11) , a(nm)
  445.       character blank   , ang(9) , out(121) , point
  446.  
  447.     1 format("1",20x,"CHART  ",i3,/)
  448.     2 format(" ",f8.3,3x,111a1)
  449.     3 format(" ")
  450.     7 format(" ",11x,".",10(9x,a1))
  451.     8 format(" ",5x,11f10.4)
  452.  
  453. * Initialize Parameters
  454.  
  455.       blank = " "
  456.       point = "."
  457.       do 5 i = 1 , 9
  458.         ang(i) = cnvi2c(i)
  459.     5 continue
  460.       nchars = nchar/10
  461.       nchars = 10*nchars
  462.       chars  = nchars
  463.       fy     = chars/10.0
  464.       ny     = fy + 1.0
  465.       nt     = nchars
  466.       nll    = nl
  467.  
  468. * Perform sort or transposition as necessary
  469.  
  470.       if (ns) 18 , 18 , 10
  471.    10 do 14 i = 1 , n
  472.         do 13 j = i , n
  473.           if ( a(i) - a(j)) 13 , 13 , 11
  474.    11     l = i - n
  475.           ll = j - n
  476.           do 12 k = 1 , m
  477.             l     = l + n
  478.             ll    = ll + n
  479.             f     = a(l)
  480.             a(l)  = a(ll)
  481.             a(ll) = f
  482.    12     continue
  483.    13   continue
  484.    14 continue
  485.       go to 18
  486.  
  487. * Test number of lines
  488.  
  489.    18 if (nll) 20 , 19 , 20
  490.    19 nll = 50
  491.    20 write(lfn,1) n0
  492.  
  493. * Find scales for base and cross variables
  494.  
  495.       xscal = (a(n)-a(1))/(float(nll-1))
  496.       m1 = n + 1
  497.       m2 = m*n
  498.       ymin = a(m1)
  499.       ymax = ymin
  500.       do 40 j = m1 , m2
  501.         if (a(j) - ymin) 26 , 28 , 28
  502.    26     ymin = a(j)
  503.    28   if (a(j) - ymax) 40 , 40 , 30
  504.    30     ymax = a(j)
  505.    40 continue
  506.       yscal  = (ymax-ymin)/chars
  507.       min    = ymin/(yscal*fy)
  508.       ynomin = fy*yscal*float(min)
  509.       if ( ymin .lt. ynomin )  ynomin = ynomin - fy*yscal
  510.       ymin   = ynomin
  511.       ymax   = ymin + chars*yscal
  512.  
  513. * Find base variable print position
  514.  
  515.       xb   = a(1)
  516.       half = xscal/2.0
  517.       l    = 1
  518.       my   = m - 1
  519.       do 80 i = 1 , nll
  520.         f = i - 1
  521.         xpr = xb + f*xscal
  522.    48   if ( a(l) .ge. xpr-half ) go to 49
  523. *         base variable too small
  524.           l = l + 1
  525.           go to 48
  526.    49   if ( a(l) .gt. xpr+half ) go to 70
  527.  
  528. * Find cross-variable
  529.  
  530.         do 55 ix = 1 , nt
  531.    55   out(ix) = blank
  532.         do 60 j = 1 , my
  533.           ll = l + j*n
  534.           jp = ((a(ll)-ymin)/yscal) + 1.5
  535.           out(jp) = ang(j)
  536.    60   continue
  537.  
  538. * Print line and clear or skip
  539.       l    = 1
  540.         write(lfn,2) xpr , ( out(iz),iz=1,nt)
  541.         l = l + 1
  542.         go to 80
  543.    70   write(lfn,3)
  544.    80 continue
  545.  
  546. * Print cross-variable numbers
  547.  
  548.       write(lfn,7) (point,ip=1,ny)
  549.       ypr(1) = ymin
  550.       do 90 kn = 1 , ny-2
  551.    90 ypr(kn+1) = ypr(kn) + yscal*chars/(ny-1)
  552.       ypr(ny) = ymax
  553.       write(lfn,8) (ypr(ip),ip=1,ny)
  554.  
  555.       return
  556.       end
  557.